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ABSTRACT This general talk is geared toward non-specialists 
and addresses three questions: (1) What types of cosmological sim- 
ulations are done? (2) What are the relevant problems to address 
and what results can one expect? (3) How can one tell if the re- 
sults are reasonable? 

1. THE SCALE OF THE PROBLEM 

The universe is huge. No, it's more like gigantically humongous. Or 
maybe it is better described as hyper-super-mega-incomparably be- 
yond any dimension we could ever imagine. And really, it's bigger than 
that, too. I sure wouldn't want to have to paint it. 

In cosmology, researchers have the audacity to attempt numerical 
simulations of the universe. Of course, one can only address a limited 
range of scales, usually referred to as the dynamic range of a simula- 
tion. The choice to focus upon certain objects or processes fixes the 
central scales in the numerical model. The push for larger models is to 
provide a wider dynamic range, perhaps extending orders of magnitude 
above and below the central scale. The goal is to fit as much as possi- 
ble of the relevant scales into a finite computational domain. 

The extent of the desired domain poses quite a challenge. From su- 
perclusters and voids and walls to clusters and groups of galaxies and 
down to giant through dwarf galaxies, cosmology attempts to explain 
structure formation on scales from kpc to hundreds of Mpc. The time 
scales involved span from Myr dynamical times to the 10 Gyr age of 
the universe. The range in mass starts at the 10^ M© scale necessary to 
resolve galaxies and extends to the 10"^^ M© scale comprising a repre- 
sentative sample of the universe. These dynamic ranges, 10^ in length, 
10"^ in time, and 10^ in mass, are the ideals and no simulation has yet 
to encompass them all. 

Similarly, the physics covers a wide range of processes. In simplest 
form, simulations model gravity via Newton's law. A coUisionless fluid 
(a.k.a. dark matter) and perhaps a cosmological constant or topological 



defects are evolved within an expanding universe. To follow the lumi- 
nous matter, Euler's equations of hydrodynamics governs the evolution 
of a coUisional baryonic fluid (a.k.a. gas). Gas dynamics is more com- 
plex and requires the inclusion of artiflcial viscosity to produce shocks 
and radiative heating and cooling. Further, the gas will eventually col- 
lapse and form stars, converting to a coUisionless stellar fluid. 

Difficulties arise when physical processes on one scale can affect 
the state of things on an extremely disparate scales, such as the radi- 
ation of a AU scale quasar ionizing gas over Mpc scale volumes. Chief 
among these effects are star formation, quasar formation, and their as- 
sociated energy, radiation, and metallicity feedbacks. As the relevant 
scales are beyond hope of inclusion in the near term, these processes 
can be included only via analytical or heuristic models. 

Considering the scope of the problem, cosmological simulations 
are one of the largest computing challenges. Cosmological codes have 
been adapted to nearly every serial, vector, or parallel supercomputer 
architecture and are an important application for testing new machines 
and driving further development. Several High Performance Computing 
and Communication grants have been awarded to groups concentrating 
on cosmological simulations, and it promises to be a forefront research 
problem well into the future. At least until the advent of hyper-super- 
mega-incomparable computers. 

2. DARK MATTER SIMULATIONS 

The basic simulation scheme evolves a coUisionless system via gravity. 
Although other methods are in use, the most common one models the 
mass as a collection of dark matter particles, also called N-body simula- 
tions. Several methods have been developed to solve Poisson's equation 
for the system; primarily by direct summation of forces between parti- 
cles (PP, Tree), by calculation of forces on a grid (PM), or by a com- 
bination of these techniques (P^M, Tree-PM). The ideas to remember 
are that grid calculation is quick, but has limited resolution, while the 
direct summation and combined methods can have arbitrarily high res- 
olution (in theory, though not in practice) at the cost of much, much 
more computational time. 

The simulation method is essentially an initial value problem. One 
chooses a cosmological model to provide a guess of the conditions in 
the early universe. A realization of these conditions for a specified vol- 
ume is then evolved forward in time. Comparing the structures in the 
simulation to observed structures provides tests of the model. Given a 
plethora of simulations and observations, of which some of each group 



are probably in error, no model fits everything. The favored model 
ends up as the one that is least inconsistent with current observations. 

Dark matter simulations are excellent for examining large scale dif- 
ferences between models. Comparisons of simulations to redshift survey 
results shifted the focus away from the top-down theory of structure 
formation, as exemplified by the hot dark matter model, and toward 
the bottom-up, or hierarchical, theory of structure formation, typified 
by the cold dark matter model. Matter distribution and correlation 
analysis provide a strong differentiation even in the mildly non-linear 
regime. 

However, in the strongly non-linear regime populated by galax- 
ies, the predictions of dark matter simulations are fuzzy. Figure 1 (all 
figures are adapted from those in Summers 1993 and Evrard, Sum- 
mers, & Davis 1994) shows a zoom into a small cluster of galaxies in 
a dark matter simulation. On the cluster scale, the dark matter halo, 
as it is called, is quite smooth with little apparent substructure. The 
one galaxy halo that can be picked out becomes quite bland upon en- 
largement. A coUisionless fluid erases substructure efficiently within 
collapsed objects. The difficulty of identifying galaxies within clusters 
as well as the assumptions one must make about how galaxies are dis- 
tributed relative to the dark matter create large error bars on the small 
scale predictions of dark matter simulations. 

That is not to say that dark matter simulations are no longer use- 
ful. If the universe is indeed dominated by coUisionless dark matter, 
then observations that are directly related to the mass are prime can- 
didates for testing models with these simulations. Such observations 
include the mass function of clusters, cluster density profiles, and gravi- 
tational lensing. Because the physics of gravity is computationally sim- 
pler and cheaper than hydrodynamics, dark matter simulations are the 
preferred method of studying these problems. 

While some of the numerics of these simulations are not germane 
to the non-specialist, a few simple checks of resolution can be made by 
anyone. The largest length scale that can appropriately be addressed in 
a simulation of box size L is roughly 0.25 L. Larger scales, and some- 
times even this scale, will be missing power from long wavelengths that 
are not included. The length resolution scale for grid based simula- 
tions is nominally given as one grid cell size (i.e., L/M for grid 
cells), but the true resolution is 2 - 3 times larger. For the particle di- 
rect summation techniques, the nominal resolution scale is the gravita- 
tional softening length, the scale over which the gravitational forces are 
smoothed to avoid point mass interactions. Typically, for particles. 



Figure 1 The dynamic range in dark matter simulations. Shown 
here are successive enlargements of a region containing a small clus- 
ter of galaxies. The top panels are slices 700 kpc thick and show 
only one fourth of the particles. Bottom panels are cubical and 
show all particles. 

the softening is set near 1/16 L/N. The true resolution depends on the 
mathematical form of the force softening, but is invariably 2-3 times 
larger than nominal. The current maximum nominal dynamic ranges 
are about 500 for grid codes and 4000 for the high resolution particle 
codes. 

Linear dynamic range is not the end of the story. The choice of 
timestep is important because one does not want particles to pass by 



each other too quickly to interact. The standard choice is to keep the 
maximum distance a particle will travel down to less than a third of a 
nominal resolution length. In equation form, it reads 



3 Vjjiax 

with At as the timestep, Ax the nominal resolution scale, and v^ax 
the maximum particle velocity. For a 10"^" year evolution with 10 kpc 
resolution and simulating 2000 km/s clusters, one needs a timestep 
At < 1.2 Myr or about 8300 timesteps (assuming an extra factor of 
to cover the tail of the velocity distribution). The dynamic range 
in mass has a nominal resolution for particle simulations of simply one 
particle mass while true resolution is of order ten particles. In grid sim- 
ulations which use particles, note that the nominal mass resolution is 
the mass within an entire grid cell, and cannot be subdivided into par- 
ticles within a cell. Any simulations which deviate strongly from the 
above ideas should justify their choices. 



3. EULERIAN HYDRODYNAMICS 

To include a baryonic component, one may model it as gaseous fluid on 
a flxed (Eulerian) grid. These codes (PPM, TVD), also called flnite dif- 
ference methods or computational fluid dynamics, solve the equations 
of hydrodynamics on the grid and couple the system to a gravity solver. 
The standard choice of gravity solver is an N-body code that also com- 
putes on a grid (PM), thereby providing similar resolution scales for 
both processes. 

The best applications of Eulerian codes are examinations of the 
gaseous component of the universe. Such applications include studies of 
the X-ray gas in clusters, Lyman a clouds, the inter-galactic medium, 
and the relative distributions of gas and dark matter. Work on X-ray 
clusters has determined that the COBE normalized cold dark matter 
model produces too many and too luminous clusters, while both the 
mixed dark matter model (some cold, some hot) and a cold dark mat- 
ter model with a cosmological constant can flt the observations ade- 
quately (e.g., Bryan et al. 1994; Kang et al. 1994). Other research has 
suggested that Lyman a clouds are not the individual collapsed objects 
their moniker would imply, but instead may be gas in an earlier stage 
of collapse, that within caustics (e.g., Cen et al. 1994). These simula- 
tions should continue to provide new ideas on gas from voids to moder- 
ate overdensities. 



And the future holds considerable promise. Computational fluid 
dynamics is a large and well sponsored fleld, having military, aerospace, 
and automotive applications. Cosmology has tapped only a portion of 
these resources. Most promising are the techniques which adapt the 
grid to the distribution of the fluid by adding more cells in high density 
regions. These techniques increase the spatial dynamic range as the 
simulation progresses and will allow one to follow the collapse process 
considerably further. 

Dynamic range considerations for hydrodynamic simulations have 
only a few differences compared to those discussed previously. Grid 
based methods have the same nominal / true length resolution of 1 / 
2-3 grid cells with mass resolution determined by the grid scale. Par- 
ticle based methods (discussed in §4 below) have their resolution set 
by the requirement that the smoothing procedure should encompass 
a scale of about 30 particles. Thus, true mass resolution is around 30 
particle masses and the spatial resolution scale is variable - larger in 
low density regions and smaller in high density regions. Nominal length 
resolution is given by the minimum smoothing length (again, true is 
2-3 times larger), which is often flxed to be near the gravitational soft- 
ening scale. The current maximum nominal dynamic ranges are about 
500 for grid codes (limited by RAM size) and 2000 for particle codes 
(limited by CPU cycles). 

The other major watchpoint for hydro codes is the physics which is 
not included in the model. "What is missing?" is an important ques- 
tion to ask of the results. Will ionizing radiation from quasars slow 
or prevent galaxy collapse (Babul & White 1991)? Does line emission 
from metals contribute to the X-ray flux in the relevant bandpass (Cen 
et al. 1995)? Is the energy injected from supernovae a signiflcant com- 
ponent in heating an intra-cluster medium (Lowenstein & Mushotzky 
1995)? These and other questions will continue to shape and reflne 
interpretation of simulation results. The bottom line is that, due to 
physics not included, some measures are meant to be considered as un- 
derestimates or as upper limits. 

4. LAGRANGIAN HYDRODYNAMICS 

The process of increasing resolution in collapsing areas happens nat- 
urally when one switches to a LaGrangian formulation (i.e., the cells 
move with the fluid). Three dimensional LaGrangian grid systems can 
produce about a factor of ten increase in range, but then break down 
when the grid becomes strongly distorted. The more successful ap- 
proach in cosmology has been the technique of smoothed particle hy- 



Figure 2 The dynamic range in SPH simulations. The same sim- 
ulation and the same panels are used as in Figure 1, but now the 
baryonic gas particles are displayed. 

drodynamics (SPH), in which gas particles act both as tracers of the 
fluid and as calculation points for the hydrodynamics. SPH uses ker- 
nel smoothing with variable smoothing lengths to cover a wide range 
of scales. As a particle based method, SPH joins cleanly with the high 
resolution particle algorithms for gravity (P^MSPH, TreeSPH). 

A direct comparison to the dark matter models is provided in Fig- 
ure 2. On the larger scales, dark matter and gas distributions are simi- 
lar because gravity is the dominant force and the collapse is essentially 
pressureless until a shock is reached. Within the collapsed regions. 
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Figure 3 An SPH disk galaxy. All baryon particles with tempera- 
ture less than 3 X 10"^ K are plotted within a 60 kpc cube. 





Group 11 55 +2526 Simulated Group 



Figure 4 The left panel shows a scanned finding chart from the 
Palomar survey for the group at 11 55, +25 26. The right panel 
shows a group from an SPH simulation. Linear scales are similar, 
roughly 500 kpc across. 



however, the baryons have cooled and condensed to galaxy scales, and 
these simulated galaxies populate the cluster. Many of these objects 
show rotationally supported disk structures (Figure 3). Groups of these 
simulated galaxies can show a striking resemblance to observed galaxy 
groups (Figure 4). 

SPH simulations are well suited to studies of highly collapsed ob- 



jects. Galaxy formation studies can delve into the formation processes, 
formation epochs, merging histories, (rough) morphologies, and other 
characteristics as a function of redshift (e.g., Katz, Hernquist, & Wein- 
berg 1993; Evrard et al. 1994). One can also follow the aggregation of 
larger structures, gaining insight into the formation processes, struc- 
ture, dynamics, and gas content of groups and clusters (e.g., Pildis, 
Bregman, & Evrard 1995; Navarro, Frenk, & White 1995; Frenk et 
al. 1995). When sufficient dynamic range is available, one can directly 
compare the statistics of simulated galaxies against large scale structure 
measures from redshift surveys. As a highlight of my recent work, the 
next section describes a premier example of SPH simulations. 

5. AN IMAX FILM PREVIEW 

In May 1996, the film Cosmic Voyage will debut at the five story 
tall IMAX theatre in the National Air and Space Museum. This film 
will take the viewer on a journey through the grand sweep of space and 
time: from the realms of sub-atomic particles to the largest scales of 
the known universe, from the present day back to the earliest seconds 
after the Big Bang. Members of the Grand Challenge Cosmology Con- 
sortium are providing the scientific basis for a four minute sequence in 
Cosmic Voyage that depicts the formation of structure in the universe. 
Beginning shortly after the Big Bang, the sequence follows the expan- 
sion of the universe, the gravitational collapse of structure, the forma- 
tion of galaxies, and a collision of of two grand design spiral galaxies. 
This sequence is the only part of the film that is direct scientific visu- 
alization: based on data from scientific simulations rather than artist's 
conceptions. 

The three parts of the sequence are being computed separately. 
Greg Bryan (Illinois) and FJS are working with the computer graphics 
experts at PIXAR to create a visualization of the early universe based 
on standard power spectrum formulae. The challenge is to coherently 
visualize and animate the tiny fluctuations in the density fleld within 
a three dimensional, expanding, and evolving universe. In the second 
part, FJS is creating a simulation of galaxy formation with unprece- 
dented dynamic range. A dramatic zoom shot will cover scales from the 
fllaments that connect clusters of galaxies down to the collapse of an 
individual galaxy. The third section features Chris Mihos' (Santa Cruz) 
simulation of two galaxies colliding and merging with intricate tidal 
streams of stars and e well as violent bursts of star formation. 

Depth scale will be added by continuing a portion of FJS' simulation as 
a background to the merger. 



Computing for IMAX resolution pushed the simulation art to the 
limits. The galaxy formation simulation (FJS) is a four million particle 
SPH run that required several CPU weeks on a dedicated 8 processor 
SGI Power Challenge. It produced 110 GB of data. The merger sim- 
ulation also used SPH techniques to follow a quarter million particles. 
750 CPU hours of computation on a Cray C90 produced over 65 GB of 
data. Both simulations are the largest of their kind performed to date. 
Visualizations are being created by a team at NCSA using proprietary 
software written at PIXAR. Today, I can show a short test video of 
the galaxy formation sequence which conveys a sense of the enormous 
amount of structure in the simulations and the stunning images that 
wiU be in the IMAX film. 

[Insert mondo-cool video clip here.] 
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